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ABSTRACT 


We  are  investigating  the  crustal  and  upper  mantle  velocity  structure  of  the  region  expressing  the  continental 
collision  between  the  Arabian  and  Eurasian  plates  using  a  joint  inversion  of  body  wave  arrival  times  and  satellite 
gravity.  The  body  wave  data  set  is  derived  from  previous  and  on-going  work  on  location  calibration  and  includes  a 
large  (-1700  events)  subset  of  events  that  qualify  as  GT590.  The  associated  arrival  time  data  sets  for  these  events 
include  many  readings  of  direct  crustal  P  and  S  phases,  as  well  as  regional  (Pn  and  Sn)  and  teleseismic  phases.  The 
data  set  has  been  carefully  groomed  to  identify  and  remove  outlier  readings  and  empirical  reading  errors  are 
estimated  for  most  arrivals  from  a  multiple  event  relocation  analysis.  The  gravity  data  provides  valuable 
supplemental  information  on  the  lateral  distribution  of  velocity  variations  that  helps  to  compensate  for  incomplete 
raypath  coverage  in  the  body  waveform  data  set,  especially  when  bandpass  filtering  is  used  to  constrain  the  gravity 
signal  to  depth  ranges  that  are  consistent  with  the  body  wave  data.  We  will  carry  out  experiments  with  smaller, 
higher-quality  GT  data  sets  as  well  as  full  data  sets  that  include  all  available  arrival  time  data  in  the  region,  and  try 
to  determine  if  the  increase  in  precision  and  accuracy  of  GT-derived  data  adequately  compensates  for  reduced 
numbers  of  data  and  raypaths.  We  will  compare  our  results  on  lateral  heterogeneity  of  the  crust  and  upper  mantle 
velocities  in  this  region  with  other  studies  and  offer  our  results  for  integration  into  Lab-based  research  efforts  to 
construct  improved  3-D  velocity  models  to  support  monitoring  efforts  in  the  region. 
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OBJECTIVES 

Our  objective  in  this  project  is  to  gain  increased  understanding  of  the  lateral  variations  of  velocity  structure  in  the 
crust  and  upper  mantle  of  the  region  expressing  the  continental  collision  between  the  Arabian  and  Eurasian  plates. 
Our  strategy  is  to  combine  different  types  of  geophysical  data  sets,  in  this  case,  body  wave  travel  time  data  and 
gravity  data,  in  a  formal  joint  inversion  for  velocity  structure.  Preliminary  work  has  shown  that  body  wave  data  and 
gravity  data  can  be  combined  to  yield  improved  resolution  of  the  variability  of  crustal  and  upper  mantle  velocities.  A 
key  element  to  this  success  is  careful  filtering  of  the  gravity  data  so  that  the  signal  source  region  corresponds  to  that 
which  is  sampled  by  the  seismic  data. 

Our  choice  of  region  for  this  study  is  driven  in  part  by  the  availability  of  an  exceptional  data  set  of  body  wave  travel 
time  data  that  results  from  a  series  of  research  projects  on  calibrated  locations  in  the  region.  No  other  area  of 
monitoring  interest  is  so  well  sampled  by  ground  truth  events.  This  project  is  one  effort  to  make  maximum  use  of 
this  unique  data  set. 

RESEARCH  ACCOMPLISHED 

This  is  a  new  project  and  efforts  so  far  have  focused  on  assembling  the  computer  codes  and  data  sets. 

Body  Wave  Data 

The  data  set  of  body  wave  arrival  times  that  will  be  used  in  this  study  will  be  derived  from  the  catalog  of 
earthquakes  in  the  study  region  that  has  been  steadily  developed  and  groomed  by  Bergman  in  collaboration  with 
E.R.  Engdahl  over  the  last  10  years  (e.g.,  Bergman  et  al.,  2009).  The  status  of  this  catalog  as  of  September  2009  is 
described  here. 
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Figure  1.  Events  in  our  catalog  that  satisfy  the  condition  of  having  secondary  azimuth  gap  <  180°,  calculated 
over  the  entire  range  of  epicentral  distances.  Secondary  azimuth  gap  is  the  largest  azimuth  gap  that 
is  filled  by  a  single  station.  Focal  depths  are  color-coded. 

•  The  catalog  contains  28,267  earthquakes  in  the  study  region,  dating  from  1923-2008.  Most  of  the  larger  and  better- 
recorded  events  have  been  reviewed  by  Engdahl  (in  single-event  location  mode,  except  for  the  most  recent  (post- 
2006)  events,  but  new  phase  readings  have  been  added  to  many  of  them. 
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The  catalog  is  substantially  complete  through  2006,  the  latest  time  for  which  we  have  ISC  data.  Data  for 
events  in  2007  and  2008  are  mainly  from  the  USGS  Earthquake  Data  Reports  (EDR)  and  local  networks 
(through  May  2008).  Some  data  from  local  stations  is  reported  by  the  EDR  in  the  second  half  of  2008, 
when  NEIC  analysts  have  manually  extracted  the  data  from  seismic  network  websites  for  events  of  interest. 

Depths:  Where  free  depth  solutions  are  not  possible  and  there  is  insufficient  depth  phase  data  (or  waveform 
studies)  to  constrain  a  fixed  depth  location,  focal  depths  are  set  at  regionalized  default  depths  based  on  the 
analysis  of  Engdahl  et  al.  (2006).  Depths  typically  range  from  about  6-20  km. 

Magnitudes:  Events  with  reported  magnitude  2.5  or  greater  are  retained  for  the  catalog,  but  some  events 
with  declared  magnitude  less  than  2.5  are  included  because  they  are  well  recorded.  A  number  of  events  in 
the  ISS  period  (1923-1963)  have  no  reported  magnitude  but  are  known  to  be  larger  events.  In  the  modern 
period,  some  events  have  no  reported  magnitude  but  have  many  phase  readings  and  are  probably  larger 
than  2.5. 

9657  events  satisfy  a  criterion  of  having  secondary  azimuth  gap  <  180°.  Events  that  don't  meet  this  criterion 
have  less  accurate  locations  and  are  probably  not  suitable  for  tomography. 

We  will  carry  out  the  joint  inversions  using  the  data  set  of -10,000  events  which  have  a  secondary  azimuth  gap  of 
less  than  180°  (Figure  1).  The  locations  and  origin  times  of  these  events  are  not  calibrated  but  their  locations  should 
have  less  bias  than  average  and  many  of  them  have  been  reviewed.  This  would  represent  a  data  set  that  is  fairly 
typical  of  those  used  in  most  tomographic  studies.  The  density  of  raypath  coverage  for  the  phase  Pn  in  an  earlier 
version  of  the  catalog  is  shown  in  Figure  2.  Coverage  in  the  southeastern  portion  of  the  region  has  been  poor 
because  of  the  scarcity  of  seismic  stations  in  that  area,  as  well  as  in  neighboring  countries.  Several  stations  have 
been  installed  in  the  region  in  recent  years  and  we  hope  for  improvements  in  data  coverage  during  the  course  of  this 
project. 


0  100  500  1000  2000  4000  8000  10000 

Figure  2.  Ray  density  map  for  Pn  in  the  study  region,  using  an  earlier  version  of  the  regional  catalog. 

We  will  also  carry  out  joint  inversions  with  a  data  set  of  much  higher  quality,  consisting  of  the  events  that  are 
members  of  calibrated  earthquake  clusters  in  the  region  (Figure  3). 
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Figure  3.  Location  map  of  35  clusters  of  earthquakes  in  the  region  that  have  been  calibrated  for  absolute 

location  and  contain  events  that  qualify  as  GT59o  or  better  (~  1700  events).  Solid  red  circles  are  all 
earthquakes  in  a  cluster,  regardless  of  calibration  level,  but  almost  all  events  are  GT1090  or  better. 

The  clusters  contain  1834  events,  of  which  1699  qualify  as  GT590.  Calibration  is  done  through  a  multiple  event 
relocation  analysis  and  the  use  of  near-source  data  of  various  types.  Origin  times  are  also  calibrated,  although  they 
are  coupled  to  focal  depth,  which  is  less  well  resolved  than  other  source  parameters.  Through  the  calibration 
analysis,  the  arrival  time  data  sets  of  these  calibrated  events  undergo  a  rigorous  and  thorough  grooming  process  that 
identifies  outliers  and  leads  to  estimates  of  empirical  reading  error  for  most  phase  arrivals. 

it  is  very  unusual  to  have  such  a  high  quality  data  set  that  provides  coverage  over  such  a  large  region.  Even  so,  the 
raypath  coverage,  especially  for  direct  crustal  phases  (Pg,  Sg),  from  the  calibrated  clusters  leaves  many  parts  of  the 
region  poorly-sampled  or  unsampled.  We  plan  to  carry  out  joint  inversions  in  more  limited  regions  in  which 
coverage  by  the  calibrated  raypaths  is  adequate.  These  results  can  be  compared  to  the  results  obtained  from  the 
larger,  uncalibrated  data  set,  to  provide  a  validation  process.  Validation  will  also  be  performed  by  calculating 
residuals  of  the  GT  dataset  through  models  derived  from  joint  inversion  of  gravity  and  uncalibrated  body  wave  data. 

We  will  also  explore  the  utility  of  performing  multiple  event  relocation  analysis  on  clusters  of  earthquakes,  even 
when  there  is  insufficient  data  to  calibrate  the  location  of  the  cluster.  This  is  likely  to  improve  resolution  in  the  joint 
inversion  analysis  in  several  ways: 

•  Relative  locations  of  all  events  in  a  cluster  will  be  significantly  improved. 

•  Events  with  poorly-constrained  locations  (even  though  they  may  pass  the  requirement  for  secondary  azimuth  gap) 
can  be  identified  and  removed  from  the  data  set. 

•  The  relocation  analysis  is  very  effective  at  identifying  outlier  readings,  which  can  be  flagged  and  removed  from 
the  data  set  for  inversion. 

There  will  still  be  unknown  bias  in  the  location  of  the  cluster,  but  the  bias  will  be  identical  for  all  events  in  the 
cluster.  Therefore  the  traditional  “event”  terms  in  tomography  can  be  replaced  by  a  much  smaller  number  of 
“cluster”  terms. 
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Even  with  these  efforts  to  provide  raypath  coverage  with  well-constrained  sources,  there  will  be  regions  that  have 
little  or  no  coverage  for  crustal  phases,  thus  limiting  what  can  be  gained  from  the  joint  inversion  concerning  crustal 
velocity  structure,  it  is  our  expectation  that  supplementing  the  body  waveform  data  with  gravity  data  will  provide 
additional  resolution  in  these  regions. 

An  Example:  Avaj 


To  illustrate  the  application  of  our  multiple  event  relocation  methodology  for  calibrated  locations  we  present  more 
detailed  information  about  one  calibrated  cluster,  Avaj  (see  Figure  3  for  location).  This  cluster  was  based  on  a  large 
(Mw  6.4),  damaging  earthquake  on  June  22,  2002  and  it’s  aftershocks.  The  cluster  of  89  events  includes  5  events 
prior  to  the  mainshock  and  aftershocks  as  recent  as  March  2008.  In  early  studies  of  this  cluster  (Walker  et  al.,  2005) 
we  used  indirect  calibration  with  one  of  the  aftershocks  that  was  reasonably  well  located  by  a  temporary  deployment 
of  seismometers,  but  this  event  was  not  very  well  linked  to  the  rest  of  the  cluster  by  readings  at  common  stations  and 
the  calibration  was  not  very  robust.  More  recently  we  obtained  readings  from  the  aftershock  survey  for  other  events 
and  more  importantly,  readings  from  more  recent  events  recorded  at  new  permanent  seismic  stations  that  provided 
enough  azimuthal  coverage  at  local  distances  to  permit  a  very  robust  calibration  using  the  direct  method  (Figure  4). 
88  events  in  the  cluster  qualify  as  GT590  or  better. 


Figure  4.  Calibrated  locations  for  the  Avaj  cluster.  Confidence  ellipses  are  shown  for  absolute  location  at 
90%  confidence  level.  Vectors  show  the  change  in  location  from  starting  locations  (single  event 
catalog  locations).  A  circle  of  5  km  radius  shown  for  reference. 

The  calibration  was  done  with  direct  P  and  S  readings  out  to  a  distance  1.7°,  using  a  single  layer  crustal  model  that 
was  adjusted  in  thickness  and  velocity  to  match  the  slopes  of  the  main  phases  and  the  Pg-Pn  crossover  distance 
(Figure  5).  It  is  evident  that  correct  phase  identification  becomes  problematic  beyond  the  cross-over  distance.  We 
investigate  the  local  distance  readings  carefully  for  every  cluster  that  will  be  calibrated  using  the  direct  method,  to 
ensure  that  the  assumed  velocity  model  adequately  accounts  for  the  arrival  time  data  and  phase  identifications,  and 
also  to  determine  where  the  assumption  of  a  1-D  model  begins  to  break  down.  When  there  is  data  at  very  short 
range,  such  as  the  S-P  readings  in  Figure  5,  the  degree  of  flattening  of  the  observed  travel  times  provides  strong 
constraint  on  focal  depth.  In  this  case  there  is  little  flattening  because  the  earthquakes  are  quite  shallow.  An  average 
assumed  depth  of  8  km  fits  the  data  well. 
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Figure  5.  Travel-time  vs.  distance  plot  for  local  distances  for  the  Avaj  cluster.  Pg  and  Sg  are  shown  as  red 

X’s.  Pn  and  Sn  arrivals  are  shown  in  green.  Theoretical  travel  times  from  the  assumed  single  layer 
crustal  model  are  shown.  Circles  are  other  phases.  In  particular  the  ones  at  distances  less  than  0.8° 
are  S-P  times  from  nearby  accelerometer  stations  with  uncalibrated  timing.  The  S-P  times  have 
been  added  to  the  corresponding  theoretical  P  time  for  proper  comparison  with  the  theoretical  S 
arrival  time. 
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Figure  6.  Travel-time  vs.  distance  plot  for  regional  distances  for  the  Avaj  cluster.  Pg  and  Sg  in  red,  Pn  and  Sn 
in  green,  P  in  black.  Theoretical  travel  times  at  less  than  15°  are  from  a  single  layer  crustal  model 
with  crustal  thickness  (47  km)  adjusted  to  lit  the  Pn  data.  Teleseismic  travel  times  are  from  akl35. 
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The  variability  in  arrival  times  of  Pn  are  shown  in  Figure  6.  The  variability  would  be  even  greater  with  the  original 
data  set.  This  data  set  has  been  carefully  edited  for  outlier  readings,  based  on  the  consistency  of  multiple 
observations  of  the  same  phase  at  the  same  station.  It  is  this  variability  that  we  are  primarily  interested  in  explaining 
with  a  3-D  model  determined  in  the  joint  inversion  (with  gravity  data)  described  below. 

Note  that  when  the  origin  times  are  established  by  near-source  data,  the  theoretical  travel  times  for  teleseismic  P  are 
too  early  by  several  seconds.  This  is  mainly  because  akl35  assumed  a  crustal  thickness  of  35  km.  Crustal  thickness 
in  the  study  region  ranges  from  45-60  km.  Therefore,  travel  times  are  increased  by  the  extra  path  length  in  the  crust. 
This  result  has  been  observed  with  virtually  every  calibrated  cluster  in  the  region,  and  is  one  of  the  reasons  that 
teleseismic  locations  have  systematic  bias. 

Gravity  Data 

Gravity  data  will  be  extracted  from  the  global  gravity  model  derived  from  the  Gravity  Recovery  and  Climate 
Experiment  (GRACE)  satellite  mission.  A  detailed  description  of  the  satellite  data  and  the  computation  of  the 
gravity  field  are  given  by  Tapley  et  al.  (2005).  The  new  gravity  field  is  more  accurate  than  previous  data  and  reveals 
previously  unobservable  tectonic  features.  Figure  7  shows  the  free-air  gravity  anomalies  for  the  proposed  study 
region.  Blue  colors  represent  gravity  highs  and  gravity  lows  are  represented  in  red. 
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Figure  7:  Free-air  gravity  observations  for  the  proposed  region  of  study.  These  anomalies  have  been 

extracted  from  the  global  gravity  model  derived  from  the  GRACE  satellite  mission.  Red  colors 
represent  gravity  lows,  blue  colors  indicate  gravity  highs.  Data  are  available  at 
http://www.csr.utexas.edu/grace/gravity,  last  accessed  June,  2010. 

Free-air  gravity  anomalies  contain  information  not  only  of  the  subsurface  density  but  also  topography.  While  in  flat 
regions  this  may  not  be  a  problem,  in  areas  of  high  topographic  relief  this  effect  needs  to  be  removed;  thus,  free-air 
anomalies  will  be  converted  into  Bouguer  anomalies  assuming  a  standard  density  for  crustal  rocks  of  2670  kg/m3. 

Modrak  et  al.  (2010)  have  shown  that  filtering  the  gravity  data  provides  important  balance  such  that  crustal  density 
variations  of  higher  spatial  frequency  can  be  exploited  for  imaging  shallow  structure,  whereas  mantle  contributions 


36 

Approved  for  public  release;  distribution  is  unlimited. 


2011  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


can  be  assumed  to  produce  the  longer-period  gravity  signals.  Judicious  application  of  filtering  is  an  important  aspect 
of  our  analysis  so  that  contributions  from  crust  and  mantle  are  handled  appropriately,  and  we  will  employ  the 
methods  of  Modrak  et  al.  (2010)  with  attention  to  selecting  the  optimal  spectral  band  to  employ  in  our  joint 
inversion. 

Joint  Inversion  of  Body  Wave  and  Gravity  Data 

The  results  of  Maceira  and  Ammon  (2009)  demonstrated  that  joint  inversion  of  surface  wave  and  gravity  data  offers 
a  simple  and  elegant  compromise  between  fitting  both  data  sets  and  improves  the  Vs  model.  In  turn,  this  will 
enhance  earthquake  location  and  characterization.  In  this  project  we  explore  the  power  of  a  similar  joint  inversion 
method,  combining  seismic  body  wave  travel  times  and  Bouguer  gravity  anomalies. 

Our  joint  inversion  of  seismic  body  wave  travel  times  (P  and  S)  and  Bouguer  gravity  anomalies  will  be  undertaken 
using  a  new  code  adapted  from  Maceira  and  Ammon  (2009)  and  integrated  into  the  regional  version  of  TomoDD 
(Zhang  and  Thurber,  2003). 

Our  treatment  of  the  gravity  data  is  similar  to  that  of  Maceira  and  Ammon  (2009).  The  gravity  anomaly  g  at  a  point 
(x,  y,  z)  is  due  to  density  anomalies  )  in  a  volume.  We  follow  Plouff  (1976)  to  calculate  the  3-D  gravity 
anomalies  of  a  prism  with  arbitrary  dimensions  (Figure  8). 


Figure  8:  Gravity  anomaly  calculation  for  a  3D  model,  (a)  The  gravity  anomaly  associated  with  a  rectangular 
prism  of  constant  density  can  be  computed  analytically  following  Plouff  (1976).  (b)  The  gravity 
anomaly  at  a  given  observation  point  Ag(m,  n)  (black  square)  can  be  computed  by  adding  the 
contribution  of  all  the  resultant  gravity  anomalies  of  all  prisms  in  the  shaded  area  of  the  model. 

In  Maceira  and  Ammon  (2009)  the  prisms  used  in  the  computation  of  the  gravity  anomalies  coincide  with  the  cells 
in  the  gridded  model.  Therefore  each  cell  in  the  model  is  a  column  of  prisms  with  horizontal  dimensions  of  1°  by  1° 
and  the  corresponding  layer  thickness  value  as  the  third  dimension.  The  contribution  of  each  cell  to  the  gravity  value 
at  a  given  point  is  a  decaying  function  of  the  distance.  Beyond  a  certain  distance,  the  effect  is  negligible.  Therefore, 
the  sensitivity  of  gravity  anomaly  to  the  density  (or,  through  scaling,  to  seismic  velocity)  is  only  calculated  for  cells 
within  a  prescribed  distance.  The  sensitivities  of  gravity  to  density  in  each  cell  are  calculated  by  perturbing  its  value 
by  5%  while  keeping  the  values  of  other  cells  fixed. 

We  have  modified  this  approach  by  attaching  the  density  values,  calculated  previously  in  each  lxl  degree  column  of 
blocks,  to  the  nodes  that  are  being  modeled  on  the  3D  grid.  This  augments  the  matrix  we  invert  for  our  updated 
model,  and  increases  the  size  of  the  problem  accordingly.  Solving  for  the  full  three-dimensional  system  has  required 
the  adoption  of  an  LSQR  approach,  rather  than  the  SVD  used  in  Maceira  and  Ammon’s  one-cell-at-a-time  approach. 

Density  -  Velocity  Relations 

One  of  the  difficulties  with  joint  inversions  is  to  determine  a  relationship  between  the  independent  data  sets.  In  this 
case,  we  require  constraints  between  seismic  velocities  and  density.  There  is  not  a  unique  and  universal  relationship 
applicable  to  all  types  of  lithologies  at  every  single  depth  under  all  possible  conditions  of  temperature  and  pressure. 
LANL  is  currently  testing,  in  collaboration  with  M.I.T.,  three  different  relationships  between  seismic  velocities  and 
density:  (1)  a  combination  of  two  existing  empirical  relationships;  one  more  suitable  for  sedimentary  rocks  after 
Nafe  and  Drake  (1963)  and  the  well-known  Birch’s  (1961)  law  more  appropriate  for  basement  rock;  (2)  Brocher’s 
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(2005);  and  (3)  Harkrider’s  (personal  communication).  An  advantage  of  the  Harkrider  relationship  is  that  it  allows 
independent  estimates  of  Vp  and  Vs  as  functions  of  density,  thus  eliminating  the  requirement  of  making  an  a  priori 
assumption  of  a  fixed  Poisson  ratio.  The  new  code,  JointTomoDD,  is  currently  being  tested  at  LANL  as  well  as  at 
the  Massachusetts  Institute  of  Technology  (MIT),  for  a  variety  of  tectonic  settings  such  as  the  Tarim  Basin,  China, 
the  crust  and  subduction  zone  beneath  western  Colombia,  and  a  thermally  active  region  within  Utah  in  the  central 
United  States. 

Preliminary  tests  of  the  body-wave  and  gravity  joint  inversion  suggest  that  a  key  influence  on  the  successful 
modeling  is  that  of  relative  weighting  of  the  two  datasets.  Our  work  in  this  project  will  emphasize  detailed 
characterization  of  the  tradeoffs  between  these  parameters  and  a  careful  comparison  of  the  obtained  models  with 
known  features  in  the  study  area,  although  the  ultimate  validation  will  be  testing  the  travel-time  predictions  of  the 
final  models  for  known  GT  paths. 

CONCLUSIONS  AND  RECOMMENDATIONS 


Because  this  project  is  just  starting  it  is  premature  to  formulate  conclusions  or  make  recommendations. 
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